! VSemiletov, 06.2011
! HeavysideFunc(r) = 1,   if r is outside of the surface
! HeavysideFunc(r) = -1,  if r is inside  of the surface
! HeavysideFunc(r) = 0,   if r is on the side where disc is not active
integer function HeavysideFunc(r)
use AcousticSurfaces, only: AcousticSurf, LeftDiscActive, RightDiscActive, &
                                nPoints, Contour3D, Center, x,y, Contour2D, &
                                CenterDir, Length, factor

implicit none
    real(8),intent(in) :: r(3)
    
    integer :: iPoint
    real(8) :: ScalarProd, localFactor, vec3D(3)
    real(8) :: vec2D(2), SurfVec(2), SurfNormal(2)
    real(8) :: ScalarProdCenter
    
    HeavysideFunc = 0
    
    vec3D(1:3) = r(1:3) - Center(1:3)

! from the left side of cylinder
    ScalarProd = sum(vec3D(1:3)*CenterDir(1:3))
    if(ScalarProd.lt.0)then
        if(LeftDiscActive)then
            HeavysideFunc = 1
        end if
        return
    end if
    
! from the right side of cylinder
    if(ScalarProd.gt.Length)then
        if(RightDiscActive)then
            HeavysideFunc = 1
        end if
        return
    end if
    
    localFactor = 1 + (factor-1.) * ScalarProd / Length
    vec2D(1) = sum(x(1:3)*vec3D(1:3)) / localFactor
    vec2D(2) = sum(y(1:3)*vec3D(1:3)) / localFactor
    
    do iPoint = 1,nPoints-1
        SurfVec(1:2) = Contour2D(1:2,iPoint) - Contour2D(1:2,iPoint+1)
        SurfNormal(1) = - SurfVec(2)
        SurfNormal(2) = SurfVec(1)
        
        ScalarProdCenter = sum(SurfNormal(1:2)*Contour2D(1:2,iPoint))
        ScalarProd = sum(SurfNormal(1:2)*(vec2D(1:2)-Contour2D(1:2,iPoint)))
        
        if(sign(1.0d0,ScalarProdCenter).eq.sign(1.0d0,ScalarProd))then
        ! 1. -> 1.0d0     SK: needs to be real for some compilers 
            HeavysideFunc = 1
            return
        end if
    end do

    HeavysideFunc = -1

end function HeavysideFunc
